// PhaseSpace interface
/////////////////////////////////////////////////////////////////////

#if !defined(PHASESPACE_H)
#define PHASESPACE_H


#include <iostream>
using namespace std;
#include <fstream>
#include <cmath>
#include <cstdlib>
#include "Function.h"


class PhaseSpace
{
friend ostream &operator<<(ostream &os, const PhaseSpace &X);

public : 
	    PhaseSpace(int InputN, Function *InputFunction);
	    PhaseSpace(void);
		virtual ~PhaseSpace();

	    void     Initialise(int InputN, Function *InputFunction);
	    void     RandomiseMomentum(double Temp = 1.0);
      //void     RandomiseMomentum(double Temp, double *mass);
	    void     IncrementPosition(double TimeStep);
      //void     IncrementPosition(double TimeStep, double *mass);
      double   IncrementPosition_observe(double TimeStep);
      double   Difference_observe(double *configuration);
	    void     CopyFinalPosition(double *x);
	    double   Energy();
      //double   Energy(double *mass);

		void     operator=(PhaseSpace &X);
		void	 SetPosition(double *x);

    void SafePositionPhaseSpace(ofstream &nameConf);

		virtual  void	IncrementMomentum(double TimeStep);
		virtual  double	PotentialEnergy(void);

		double	 *Position;

private:
		bool     DMA;     // DMA?
		int 	 N; 
		Function *SampleFunction;
		double   *Momentum, *Derivative;

    double   KineticEnergy();
		//double   KineticEnergy(double *mass);

    // to calculate length taken by one Nmd step

    double mLength, mIncreaseBy;    
};

//////////////////////////////////////////////////////////////////
// Phase Space : inline constructors/destructors
//////////////////////////////////////////////////////////////////


inline PhaseSpace::PhaseSpace(int InputN, Function *InputFunction)
{
Initialise(InputN,InputFunction);	
}

inline PhaseSpace::PhaseSpace(void) 
{
DMA = false;
}

//////////////////////////////////////////////////////////////////
// Phase Space : inline overloaded operators
//////////////////////////////////////////////////////////////////
// sets the position in memory
inline void   PhaseSpace::SetPosition(double *x)
{
CopyA(Position, x, N);
}


//////////////////////////////////////////////////////////////////
// Phase Space : inline functions
//////////////////////////////////////////////////////////////////
// copies the final position to memory
inline void   PhaseSpace::CopyFinalPosition(double *x)
{
CopyA(x, Position, N);
}

// randomly initialises the momentum
inline void   PhaseSpace::RandomiseMomentum(double Temp)
{
	double Scale = sqrt(Temp);

  for(int n=1;n<=N;n++)
  {
    Momentum[n] = Scale*gasdev();
  }
}

/*inline void   PhaseSpace::RandomiseMomentum(double Temp, double *mass)
{
	double Scale = sqrt(Temp);

  for(int n=1;n<=N;n++)
  {
    Momentum[n] = sqrt(mass[n])*Scale*gasdev();
    //Momentum[n] = Scale*gasdev();
  }
}*/

// calculates the potential energy
inline double PhaseSpace::PotentialEnergy(void)
{
return ( SampleFunction->Value(Position) );
}

// calculates the total energy
inline double PhaseSpace::Energy()
{
return ( PotentialEnergy() + KineticEnergy() );
}

/*inline double PhaseSpace::Energy(double *mass)
{
return ( PotentialEnergy() + KineticEnergy(mass) );
}*/



#endif // if !defined(PHASESPACE_H)
